This script analyzes and plots data for Symbiont Integration 2020 respirometry data. Plots are displayed in Plotting section. Results are provided in Analysis section along with summaries of potential implications.

Setup

Set up workspace, set options, and load required packages.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse') 
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car') 
if ("lme4" %in% rownames(installed.packages()) == 'FALSE') install.packages('lme4') 
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest') 
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales') 
if ("cowplot" %in% rownames(installed.packages()) == 'FALSE') install.packages('cowplot') 
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2') 
if ("effects" %in% rownames(installed.packages()) == 'FALSE') install.packages('effects') 
if ("emmeans" %in% rownames(installed.packages()) == 'FALSE') install.packages('emmeans') 
if ("multcomp" %in% rownames(installed.packages()) == 'FALSE') install.packages('multcomp') 

#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('lme4')
library('lmerTest')
library('scales')
library('cowplot')
library('effects')
library('emmeans')
library('multcomp')

Data visualization and manipulation

Load data from LoLinR.

PRdata<-read.csv("Mcap2020/Output/Respiration/oxygen_P_R_calc.csv") #load data

Separate project specific data.

#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]

#format columns
PRdata$dpf<-as.factor(PRdata$dpf)

#subset fused juvenile data frame for a different project, then save as separate file
fused_oxygen <- PRdata[which(PRdata$Lifestage=='Juvenile'),]
fused_oxygen <- droplevels(fused_oxygen)

Calculate a P:R ratio using gross photosynthesis.

PRdata$ratio<-abs(PRdata$GP.nmol.org.min)/abs(PRdata$R.nmol.org.min) #calculate ratio with absolute values
#remove outliers detected by values of P:R ratio data
PRdata<-PRdata%>%filter(ratio < 15)

fused_oxygen$ratio<-abs(fused_oxygen$GP.nmol.org.min)/abs(fused_oxygen$R.nmol.org.min)

Analysis and plotting of core dataset

Load sample metadata to identify stage and hpf.

metadata<-read.csv("Mcap2020/Data/Respiration/PR_metadata.csv")

PRdf<-PRdata%>%
  filter(!Lifestage=="Juvenile")

#load in sample information
PRdf$code<-paste(PRdf$Date, PRdf$Lifestage)
metadata$code<-paste(metadata$Date, metadata$Lifestage)

PRdf$hpf<-metadata$hpf[match(PRdf$code, metadata$code)]
PRdf$group<-metadata$group[match(PRdf$code, metadata$code)]
PRdf$hpf<-as.factor(PRdf$hpf)

Respiration

Plot data as a scatterplot

r_plot<-PRdf %>%
    ggplot(., aes(x = hpf, y = R.nmol.org.min)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
    #geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot

#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30

#ggsave("Mcap2020/Figures/Respiration/Respiration.png", r_plot, dpi=300, w=8.5, h=5, units="in")

Photosynthesis (Net)

Plot data as a scatterplot

p_plot<-PRdf %>%
    ggplot(., aes(x = hpf, y = P.nmol.org.min)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
    #geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
    ylab(expression(bold(paste("P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); p_plot

#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30

#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", rp_plot, dpi=300, w=8.5, h=5, units="in")

Photosynthesis (Gross)

Plot data as a scatterplot

gp_plot<-PRdf %>%
    ggplot(., aes(x = hpf, y = GP.nmol.org.min)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
    #geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position=position_jitterdodge(0.2)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
    scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
    ylab(expression(bold(paste("GP (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); gp_plot

#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30

#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", gp_plot, dpi=300, w=8.5, h=5, units="in")

P:R Ratio

Plot data as a scatterplot

pr_plot<-PRdf %>%
    ggplot(., aes(x = hpf, y = ratio)) +
    geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
    #geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
    geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) + 
    xlab("Hours Post-Fertilization") + 
    scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
    scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
    ylab(expression(bold(paste("P:R")))) +
    scale_y_continuous(limits=c(-0.5, 10), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=16),
      legend.text=element_text(size=14)
      ); pr_plot

 #EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30

#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", gp_plot, dpi=300, w=8.5, h=5, units="in")

Generate final figure

full_fig<-plot_grid(r_plot, p_plot, gp_plot, pr_plot, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")

ggsave(filename="Mcap2020/Figures/Respiration/respirometry_fig.png", plot=full_fig, dpi=500, width=26, height=6, units="in")

Analyze data

Respiration

Build linear mixed effect model and examine for Respiration.

Rmodel1<-lmer(R.nmol.org.min~group + (1|Run), data=PRdf) #run as random
anova(Rmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##         Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## group 0.024176 0.0048353     5 95.206  30.327 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Rmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: R.nmol.org.min ~ group + (1 | Run)
##    Data: PRdf
## 
## REML criterion at convergence: -551.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1454 -0.3542 -0.0011  0.4361  1.8624 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Run      (Intercept) 3.355e-05 0.005792
##  Residual             1.594e-04 0.012627
## Number of obs: 103, groups:  Run, 4
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -0.013525   0.004300  6.419042  -3.145   0.0182 *  
## groupLarvae3   0.003200   0.003993 93.407620   0.801   0.4250    
## groupLarvae5  -0.006146   0.004347 96.879347  -1.414   0.1607    
## groupLarvae6  -0.004662   0.004446 96.126246  -1.049   0.2970    
## groupRecruit1 -0.041897   0.004367 96.475039  -9.594 1.06e-15 ***
## groupRecruit2 -0.024674   0.004736 95.896340  -5.210 1.08e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.464                            
## groupLarva5 -0.538  0.459                     
## groupLarva6 -0.478  0.449  0.476              
## groupRecrt1 -0.540  0.457  0.531  0.467       
## groupRecrt2 -0.449  0.422  0.447  0.434  0.439

Check assumptions of model for residual normality and variance. Passes assumptions.

qqPlot(residuals(Rmodel1))

## [1] 49 94
leveneTest(residuals(Rmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)   
## group  5  3.4255 0.00681 **
##       97                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conduct post hoc test.

emm<-emmeans(Rmodel1, ~group)
cld(emm)
##  group     emmean      SE    df lower.CL  upper.CL .group
##  Recruit1 -0.0554 0.00416  7.70  -0.0651 -0.045773  1    
##  Recruit2 -0.0382 0.00481 12.00  -0.0487 -0.027719   2   
##  Larvae5  -0.0197 0.00416  7.70  -0.0293 -0.010022    3  
##  Larvae6  -0.0182 0.00453  9.48  -0.0283 -0.008029    3  
##  Larvae2  -0.0135 0.00441  7.74  -0.0237 -0.003305    3  
##  Larvae3  -0.0103 0.00441  7.74  -0.0205 -0.000105    3  
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

Net Photosynthesis

Build linear mixed effect model and examine for Photosynthesis

Pmodel1<-lmer(P.nmol.org.min~group + (1|Run), data=PRdf) #run nested within day
anova(Pmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## group 0.0098599 0.001972     5 94.852  6.7341 2.065e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Pmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: P.nmol.org.min ~ group + (1 | Run)
##    Data: PRdf
## 
## REML criterion at convergence: -493.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2873 -0.3392 -0.0035  0.5404  2.7571 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Run      (Intercept) 3.251e-05 0.005701
##  Residual             2.928e-04 0.017112
## Number of obs: 103, groups:  Run, 4
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)    0.012253   0.005092  8.055364   2.406  0.04255 * 
## groupLarvae3   0.002363   0.005411 93.362976   0.437  0.66341   
## groupLarvae5   0.019696   0.005831 94.930217   3.378  0.00106 **
## groupLarvae6   0.015380   0.005998 96.884244   2.564  0.01187 * 
## groupRecruit1 -0.009215   0.005849 93.503867  -1.576  0.11850   
## groupRecruit2  0.000326   0.006392 96.717807   0.051  0.95943   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.531                            
## groupLarva5 -0.598  0.464                     
## groupLarva6 -0.539  0.451  0.473              
## groupRecrt1 -0.600  0.463  0.521  0.466       
## groupRecrt2 -0.505  0.423  0.444  0.429  0.437

Check assumptions of model for residual normality and variance. Violates variance assumption, return to this.

qqPlot(residuals(Pmodel1))

## [1] 72 45
leveneTest(residuals(Pmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  5  3.2159 0.009928 **
##       97                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conduct post hoc test.

emm<-emmeans(Pmodel1, ~group)
cld(emm)
##  group     emmean      SE    df  lower.CL upper.CL .group
##  Recruit1 0.00304 0.00494 11.73 -0.007763   0.0138  1    
##  Larvae2  0.01225 0.00531  9.99  0.000419   0.0241  12   
##  Recruit2 0.01258 0.00591 18.43  0.000184   0.0250  12   
##  Larvae3  0.01462 0.00531  9.99  0.002782   0.0265  123  
##  Larvae6  0.02763 0.00548 13.93  0.015871   0.0394   23  
##  Larvae5  0.03195 0.00494 11.73  0.021148   0.0427    3  
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

Gross Photosynthesis

Build linear mixed effect model and examine for GP

GPmodel1<-lmer(GP.nmol.org.min~group + (1|Run), data=PRdf) #run nested within day
anova(GPmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##         Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## group 0.015016 0.0030032     5 95.147  6.5732 2.701e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(GPmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GP.nmol.org.min ~ group + (1 | Run)
##    Data: PRdf
## 
## REML criterion at convergence: -448.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8197 -0.4243  0.0855  0.4891  4.8386 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Run      (Intercept) 0.0001121 0.01059 
##  Residual             0.0004569 0.02138 
## Number of obs: 103, groups:  Run, 4
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    0.0256301  0.0075645  5.8288704   3.388 0.015378 *  
## groupLarvae3  -0.0008371  0.0067594 93.2900454  -0.124 0.901704    
## groupLarvae5   0.0259018  0.0073748 96.9863051   3.512 0.000677 ***
## groupLarvae6   0.0195399  0.0075337 95.8607643   2.594 0.010983 *  
## groupRecruit1  0.0329184  0.0074107 96.7726143   4.442 2.37e-05 ***
## groupRecruit2  0.0244979  0.0080232 95.6295932   3.053 0.002930 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.447                            
## groupLarva5 -0.521  0.458                     
## groupLarva6 -0.461  0.449  0.477              
## groupRecrt1 -0.523  0.456  0.533  0.467       
## groupRecrt2 -0.433  0.421  0.448  0.435  0.439

Check assumptions of model for residual normality and variance. Violates variance assumption, return to this.

qqPlot(residuals(GPmodel1))

## [1] 49 94
leveneTest(residuals(GPmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  5  4.3752 0.001236 **
##       97                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conduct post hoc test.

emm<-emmeans(GPmodel1, ~group)
cld(emm)
##  group    emmean      SE    df lower.CL upper.CL .group
##  Larvae3  0.0248 0.00772  7.24  0.00665   0.0429  1    
##  Larvae2  0.0256 0.00772  7.24  0.00749   0.0438  1    
##  Larvae6  0.0452 0.00792  8.65  0.02714   0.0632  12   
##  Recruit2 0.0501 0.00839 10.81  0.03163   0.0686   2   
##  Larvae5  0.0515 0.00731  7.04  0.03426   0.0688   2   
##  Recruit1 0.0585 0.00731  7.04  0.04127   0.0758   2   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

P:R

Build linear mixed effect model and examine for GP

PRmodel1<-lmer(ratio~group + (1|Run), data=PRdf) #run nested within day
anova(PRmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##       Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## group 36.271  7.2542     5 95.046  4.6091 0.0008271 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(PRmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ratio ~ group + (1 | Run)
##    Data: PRdf
## 
## REML criterion at convergence: 336.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5330 -0.4931 -0.0433  0.2952  6.0621 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Run      (Intercept) 0.009501 0.09747 
##  Residual             1.573887 1.25455 
## Number of obs: 103, groups:  Run, 4
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    2.51163    0.28853 21.48307   8.705 1.72e-08 ***
## groupLarvae3  -0.36950    0.39672 95.33610  -0.931  0.35401    
## groupLarvae5   0.19131    0.41028 87.81805   0.466  0.64217    
## groupLarvae6   0.00286    0.43017 95.29471   0.007  0.99471    
## groupRecruit1 -1.38807    0.41030 87.46480  -3.383  0.00107 ** 
## groupRecruit2 -1.15704    0.45965 96.02507  -2.517  0.01348 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.687                            
## groupLarva5 -0.683  0.483                     
## groupLarva6 -0.646  0.461  0.454              
## groupRecrt1 -0.683  0.483  0.480  0.454       
## groupRecrt2 -0.604  0.432  0.425  0.405  0.424

Check assumptions of model for residual normality and variance. Meets assumptions.

qqPlot(residuals(PRmodel1))

## [1] 19 11
leveneTest(residuals(PRmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  1.8893 0.1032
##       97

Conduct post hoc test.

emm<-emmeans(PRmodel1, ~group)
cld(emm)
##  group    emmean    SE   df lower.CL upper.CL .group
##  Recruit1   1.12 0.300 45.3    0.519     1.73  1    
##  Recruit2   1.35 0.371 53.6    0.611     2.10  12   
##  Larvae3    2.14 0.312 15.9    1.480     2.80  12   
##  Larvae2    2.51 0.312 15.9    1.849     3.17   2   
##  Larvae6    2.51 0.334 41.5    1.841     3.19   2   
##  Larvae5    2.70 0.300 45.3    2.099     3.31   2   
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

Generate summary tables of mean and SE for all variables

Generate summary of all respiration data.

summary<-PRdf%>%
  group_by(group)%>%
  summarise(N=length(R.nmol.org.min),
            Mean_R=mean(R.nmol.org.min), 
            SE_R=sd(R.nmol.org.min)/sqrt(length(R.nmol.org.min)),
            Mean_P=mean(P.nmol.org.min), 
            SE_P=sd(P.nmol.org.min)/sqrt(length(P.nmol.org.min)),
            Mean_GP=mean(GP.nmol.org.min), 
            SE_GP=sd(GP.nmol.org.min)/sqrt(length(GP.nmol.org.min)),
            Mean_PR=mean(ratio), 
            SE_PR=sd(ratio)/sqrt(length(ratio)))

summary%>%
  write_csv(., "Mcap2020/Output/Respiration/mean_respiration.csv")

Exploring additional datasets

Plotting: Larvae, Recruits, Juveniles

Generate plots by lifestage and days post fertilization.

Respiration

Generate dot plot.

Rplot1<-ggplot(data=PRdata, aes(x=dpf, y=R.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
  scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Days post fertilization")+
  ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "top",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Rplot1

Rplot2<-Rplot1+theme(legend.position="none")

#ggsave("Mcap2020/Figures/Respiration/mcap_R.pdf", plot=Rplot1, height=8, width=7, units = c("in"), dpi=300) #output figure

Generate mean table

meanR <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
                   N    = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
                   mean = mean(R.nmol.org.min, na.rm=TRUE),
                   sd   = sd(R.nmol.org.min, na.rm=TRUE),
                   se   = sd / sqrt(N)
); meanR
##   Lifestage dpf  N        mean          sd           se
## 1    Embryo   3 20 -0.01684416 0.008734091 0.0019530020
## 2  Juvenile  40 20 -0.04521002 0.019598294 0.0043823118
## 3    Larvae   4 20 -0.01364444 0.004434086 0.0009914917
## 4    Larvae   8 18 -0.01967560 0.005887039 0.0013875884
## 5    Larvae  10 15 -0.02028284 0.008694612 0.0022449392
## 6   Recruit   8 18 -0.05541750 0.027867296 0.0065683846
## 7   Recruit  10 12 -0.04029479 0.006264046 0.0018082742

Net Photosynthesis

Pplot1<-ggplot(data=PRdata, aes(x=dpf, y=P.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
  scale_y_continuous(limits=c(-0.1, 0.3), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Days post fertilization")+
  ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); Pplot1

Pplot2<-Pplot1+theme(legend.position="none")

#ggsave("Mcap2020/Figures/Respiration/mcap_P.pdf", plot=Pplot2, height=8, width=7, units = c("in"), dpi=300) #output figure

Generate mean table.

meanP <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
                   N    = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
                   mean = mean(P.nmol.org.min, na.rm=TRUE),
                   sd   = sd(P.nmol.org.min, na.rm=TRUE),
                   se   = sd / sqrt(N)
); meanP
##   Lifestage dpf  N        mean          sd          se
## 1    Embryo   3 20 0.011275224 0.009282584 0.002075649
## 2  Juvenile  40 20 0.079503148 0.035835414 0.008013042
## 3    Larvae   4 20 0.013637833 0.014767901 0.003302203
## 4    Larvae   8 18 0.032156372 0.011762116 0.002772357
## 5    Larvae  10 15 0.029064523 0.018198588 0.004698855
## 6   Recruit   8 18 0.002830864 0.028787115 0.006785188
## 7   Recruit  10 12 0.014010579 0.016925897 0.004886086

Gross Photosynthesis

GPplot1<-ggplot(data=PRdata, aes(x=dpf, y=GP.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
  scale_y_continuous(limits=c(-0.1, 0.3), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Days post fertilization")+
  ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); GPplot1

GPplot2<-GPplot1+theme(legend.position="none") 

#ggsave("Mcap2020/Figures/Respiration/mcap_GP.pdf", plot=GPplot2, height=8, width=7, units = c("in"), dpi=300) #output figure

Generate mean table.

meanGP <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
                   N    = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
                   mean = mean(GP.nmol.org.min, na.rm=TRUE),
                   sd   = sd(GP.nmol.org.min, na.rm=TRUE),
                   se   = sd / sqrt(N)
); meanGP
##   Lifestage dpf  N       mean         sd          se
## 1    Embryo   3 20 0.02811939 0.01377244 0.003079611
## 2  Juvenile  40 20 0.12471317 0.04749814 0.010620907
## 3    Larvae   4 20 0.02728227 0.01561791 0.003492270
## 4    Larvae   8 18 0.05183197 0.01398628 0.003296597
## 5    Larvae  10 15 0.04934737 0.02530702 0.006534244
## 6   Recruit   8 18 0.05824836 0.03865966 0.009112168
## 7   Recruit  10 12 0.05430537 0.01733436 0.005003998

Gross Photosynthesis:Respiration ratio 

PRplot1<-ggplot(data=PRdata, aes(x=dpf, y=ratio, colour=Lifestage, group=interaction(Lifestage, dpf))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
  geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-2, 6), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  ylab(expression(bold(paste("P : R ")))) +
  xlab("Days Post Fertilization") +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "top",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=22), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); PRplot1

PRplot2<-PRplot1+theme(legend.position="right") 

#ggsave("Mcap2020/Figures/Respiration/mcap_PR.pdf", plot=PRplot1, height=8, width=8, units = c("in"), dpi=300) #output figure

Generate mean table.

meanPR <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
                   N    = length(ratio[!is.na(ratio)]),
                   mean = mean(ratio, na.rm=TRUE),
                   sd   = sd(ratio, na.rm=TRUE),
                   se   = sd / sqrt(N)
); meanPR
##   Lifestage dpf  N     mean        sd        se
## 1    Embryo   3 20 2.508994 2.3952025 0.5355836
## 2  Juvenile  40 20 3.002808 1.0719306 0.2396910
## 3    Larvae   4 20 2.139497 1.0377752 0.2320536
## 4    Larvae   8 18 2.705802 0.7656908 0.1804751
## 5    Larvae  10 15 2.518476 0.5820499 0.1502846
## 6   Recruit   8 18 1.120693 0.6434859 0.1516711
## 7   Recruit  10 12 1.358576 0.4359740 0.1258549

Plot Panels

Combine plots above into a single panel.

mcap_respiration<-plot_grid(Rplot2, Pplot2, GPplot2, PRplot2, labels = c("A", "B", "C", "D"), label_size=18, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(0.75,0.75,0.75,1), align="h")

#ggsave(filename="Mcap2020/Figures/Respiration/mcap_PR_panel.pdf", plot=mcap_respiration, dpi=300, width=24, height=6, units="in")

Plotting: Fused juveniles

Polyp Specific Rates

Generate a plot for fused corals separately to examine respirometry normalized to the number polyps fused together within a colony.

Respiration

FuseR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=R.nmol.org.min, group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.1, 0.01), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("R (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseR1

Net Photosynthesis

FuseP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=P.nmol.org.min, group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.1, 0.4), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("Net P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseP1

Gross Photosynthesis

FuseGP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=GP.nmol.org.min, group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.1, 0.4), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("Gross P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseGP1

P:R ratio

FusePR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=ratio, group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-1, 8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("P : R")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FusePR1

Colony Total Rates

Respiration

Fuse_c_R1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(R.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.4, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("R (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_R1

Net Photosynthesis

Fuse_c_P1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(P.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.1, 0.8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("Net P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_P1

Gross Photosynthesis

Fuse_c_GP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(GP.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-0.1, 0.8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("Gross P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_GP1

P:R Ratio 

Fuse_c_PR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(abs(fused_oxygen$GP.nmol.org.min*Org.Number)/abs(fused_oxygen$R.nmol.org.min*Org.Number)), group=interaction(Org.Number))) +
  geom_boxplot(position=position_dodge(0.9), lwd=1)+
  geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
  geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-1, 10), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("P : R")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_PR1

Plot Panel

Combine plots above into a single panel.

mcap_fused_respiration<-plot_grid(FuseR1, FuseP1, FuseGP1, FusePR1, Fuse_c_R1, Fuse_c_P1, Fuse_c_GP1, Fuse_c_PR1, labels = c("Polyp - Level Responses", "", "", "", "Colony - Level Responses", "", "", ""), label_size = 18, ncol=4, nrow=2, rel_heights= c(1), rel_widths = c(1), align="h")

ggsave(filename="Mcap2020/Figures/Respiration/mcap_fused_panel.pdf", plot=mcap_fused_respiration, dpi=300, width=24, height=12, units="in")

Comparing colony & polyp R and P rates

Respiration

Make a summary table to displays means and standard errors for plotting.

colony_v_polypsR <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
                   N_r_polyps    = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
                   mean_r_polyps = mean(R.nmol.org.min, na.rm=TRUE),
                   sd_r_polyps   = sd(R.nmol.org.min, na.rm=TRUE),
                   se_r_polyps   = sd_r_polyps / sqrt(N_r_polyps), 
                   N_r_colony    = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
                   mean_r_colony = mean((Org.Number*R.nmol.org.min), na.rm=TRUE),
                   sd_r_colony   = sd(Org.Number*R.nmol.org.min, na.rm=TRUE),
                   se_r_colony   = sd_r_colony / sqrt(N_r_colony)
)

Plot a line plot with two axis. Left axis: polyp specific respiration. Right axis: colony total respiration.

R_c_vs_p<-ggplot(data=colony_v_polypsR, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
  geom_point(aes(y=mean_r_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
  geom_errorbar(aes(ymin=mean_r_polyps-se_r_polyps, ymax=mean_r_polyps+se_r_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
   geom_point(aes(y=mean_r_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
  geom_errorbar(aes(ymin=mean_r_colony-se_r_colony, ymax=mean_r_colony+se_r_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(
     name = "Polyp Respiration", # Features of the first axis
     sec.axis = sec_axis(~.*1, name="Colony Total Respiration"),  # Add a second axis and specify its features
     limits=c(-0.4, 0.1),
     labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
     ) + 
  xlab("Number of Fused Polyps")+
  #ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"), 
        axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));R_c_vs_p

Net Photosynthesis

Make a summary table to displays means and standard errors for plotting.

colony_v_polypsP <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
                   N_p_polyps    = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
                   mean_p_polyps = mean(P.nmol.org.min, na.rm=TRUE),
                   sd_p_polyps   = sd(P.nmol.org.min, na.rm=TRUE),
                   se_p_polyps   = sd_p_polyps / sqrt(N_p_polyps), 
                   N_p_colony    = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
                   mean_p_colony = mean((Org.Number*P.nmol.org.min), na.rm=TRUE),
                   sd_p_colony   = sd(Org.Number*P.nmol.org.min, na.rm=TRUE),
                   se_p_colony   = sd_p_colony / sqrt(N_p_colony)
)

Plot a line plot with two axis. Left axis: polyp specific net photosynthesis Right axis: colony total net photosynthesis.

P_c_vs_p<-ggplot(data=colony_v_polypsP, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
  geom_point(aes(y=mean_p_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
  geom_errorbar(aes(ymin=mean_p_polyps-se_p_polyps, ymax=mean_p_polyps+se_p_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
   geom_point(aes(y=mean_p_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
  geom_errorbar(aes(ymin=mean_p_colony-se_p_colony, ymax=mean_p_colony+se_p_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(
     name = "Polyp Net Photosynthesis", # Features of the first axis
     sec.axis = sec_axis(~.*1, name="Colony Total Net Photosynthesis"),  # Add a second axis and specify its features
     limits=c(-0.1, 0.8),
     labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
     ) + 
  xlab("Number of Fused Polyps")+
  #ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"), 
        axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));P_c_vs_p

Gross Photosynthesis

Make a summary table to displays means and standard errors for plotting.

colony_v_polypsGP <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
                   N_gp_polyps    = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
                   mean_gp_polyps = mean(GP.nmol.org.min, na.rm=TRUE),
                   sd_gp_polyps   = sd(GP.nmol.org.min, na.rm=TRUE),
                   se_gp_polyps   = sd_gp_polyps / sqrt(N_gp_polyps), 
                   N_gp_colony    = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
                   mean_gp_colony = mean((Org.Number*GP.nmol.org.min), na.rm=TRUE),
                   sd_gp_colony   = sd(Org.Number*GP.nmol.org.min, na.rm=TRUE),
                   se_gp_colony   = sd_gp_colony / sqrt(N_gp_colony)
)

Plot a line plot with two axis. Left axis: polyp specific net photosynthesis Right axis: colony total net photosynthesis.

GP_c_vs_p<-ggplot(data=colony_v_polypsGP, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
  geom_point(aes(y=mean_gp_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
  geom_errorbar(aes(ymin=mean_gp_polyps-se_gp_polyps, ymax=mean_gp_polyps+se_gp_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
   geom_point(aes(y=mean_gp_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
  geom_errorbar(aes(ymin=mean_gp_colony-se_gp_colony, ymax=mean_gp_colony+se_gp_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
  geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(
     name = "Polyp Gross Photosynthesis", # Features of the first axis
     sec.axis = sec_axis(~.*1, name="Colony Total Gross Photosynthesis"),  # Add a second axis and specify its features
     limits=c(-0.1, 0.8),
     labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
     ) + 
  xlab("Number of Fused Polyps")+
  #ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"), 
        axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));GP_c_vs_p

P:R

P:R is the same when calculated for polyp or colony rates. Display P:R means in plot.

Generate mean plot.

meanPR_fuse <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
                   N    = length(ratio[!is.na(ratio)]),
                   mean = mean(ratio, na.rm=TRUE),
                   sd   = sd(ratio, na.rm=TRUE),
                   se   = sd / sqrt(N)
)
PRplot_fuse<-ggplot(data=meanPR_fuse, aes(x=as.factor(Org.Number), y=mean, group=interaction(Org.Number))) +
  geom_point(color="darkgray", size=3, position=position_dodge(0.3)) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0, linetype="solid", position=position_dodge(0.3), size=1) +
  geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
  scale_y_continuous(limits=c(-1, 8), labels = scales::number_format(accuracy = 1, decimal.mark = '.'))+
  xlab("Number of Fused Polyps")+
  ylab(expression(bold(paste("Gross Photosynthesis: Respiration")))) +
  theme_classic() +
  theme(text = element_text(size = 18, color="black"),
        axis.text = element_text(size = 18, color="black"), 
        legend.position = "right",
        axis.title = element_text(size = 18, color="black", face="bold"), 
        legend.title=element_blank(), 
        legend.text = element_text(size=18), 
        plot.margin=unit(c(1,1,1,1), "cm"), 
        axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)), 
        axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));PRplot_fuse

Plot Panel

Combine plots above into a single panel.

mcap_fused_colony_vs_polyp<-plot_grid(R_c_vs_p, P_c_vs_p, GP_c_vs_p, PRplot_fuse, labels = c("Respiration", "Net Photosynthesis", "Gross Photosynthesis", "P:R"), label_size = 16, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,0.8), align="h")

ggsave(filename="Mcap2020/Figures/Respiration/mcap_fused_colony_vs_polyp_panel.pdf", plot=mcap_fused_colony_vs_polyp, dpi=300, width=24, height=6, units="in")

Analysis: Larvae, Recruits, Juveniles

Respiration

Build linear mixed effect model and examine for Respiration.

Rmodel1<-lmer(R.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata) #run nested within day
anova(Rmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq   Mean Sq NumDF   DenDF F value    Pr(>F)    
## Lifestage     0.0181943 0.0060648     3   6.737 32.9388 0.0002107 ***
## dpf           0.0012575 0.0006287     2   5.001  3.4148 0.1161237    
## Lifestage:dpf 0.0009344 0.0009344     1 108.357  5.0747 0.0262915 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Variance is violated, revisit this assumption.

qqPlot(residuals(Rmodel1))

## [1] 69 38
hist(residuals(Rmodel1))

leveneTest(residuals(Rmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   6  3.9085 0.001361 **
##       116                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For respiration, lifestage and lifestage:days post fertilization are significant. Respiration increased over larval development and was higher for recruits at day 8. This is the closest timepoint to metamorphosis, which may suggest that respiratory demand is high immediately after settlement and then reduces by day 10. Juveniles at 40 days post settlement had high respiration, indicating higher demand with larger size.

Net Photosynthesis

Build linear mixed effect model and examine for Net Photosynthesis.

Pmodel1<-lmer(P.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(Pmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## Lifestage     0.046409 0.0154695     3 103.74 35.1807 9.104e-16 ***
## dpf           0.002701 0.0013503     2 114.84  3.0708   0.05021 .  
## Lifestage:dpf 0.000715 0.0007147     1 113.16  1.6253   0.20496    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Variance is just barely violated, revisit this assumption.

qqPlot(residuals(Pmodel1))

## [1] 21 92
hist(residuals(Pmodel1))

leveneTest(residuals(Pmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)    
## group   6  5.0211 0.000131 ***
##       116                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For net photosynthesis, lifestage is significant. Photosynthesis appears to be higher in larvae (esp. at day 8), but reduces in early recruits (day 8 and day 10), which may suggest that symbiosis is “upset” or “reset” after metamorphosis and takes time to build to high levels seen in 40 day old juveniles. This supports our hypotheses that symbiosis integration occurs after settlement during early juvenile development.

Gross Photosynthesis

Build linear mixed effect model and examine for Gross Photosynthesis.

GPmodel1<-lmer(GP.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(GPmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## Lifestage     0.076694 0.0255645     3 111.63 36.2714 < 2.2e-16 ***
## dpf           0.007076 0.0035380     2 115.53  5.0198  0.008119 ** 
## Lifestage:dpf 0.000022 0.0000217     1 112.90  0.0308  0.860913    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance.

qqPlot(residuals(GPmodel1))

## [1] 38 69
hist(residuals(GPmodel1))

leveneTest(residuals(GPmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   6  4.3982 0.0004844 ***
##       116                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For gross photosynthesis, lifestage is significant. Gross photosynthesis does not appear to reduce as much as net photosynthesis after settlement, so this may indicate that the heavy respiratory demand is what is leading to decreased net photosynthetic output and therefore less available for growth/calcification at this life stage. Gross photosynthesis is still highest in juveniles, which supports that symbiosis is active and integrated in older juveniles.

P:R

Build linear mixed effect model and examine for P:R ratio (gross photosynthesis : respiration)

PRmodel1<-lmer(ratio~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(PRmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF   DenDF F value   Pr(>F)   
## Lifestage     39.383 13.1277     3  12.025  8.8875 0.002233 **
## dpf            2.247  1.1234     2   9.901  0.7605 0.492890   
## Lifestage:dpf  0.665  0.6654     1 109.312  0.4505 0.503519   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions met.

qqPlot(residuals(PRmodel1))

## [1] 19 11
hist(residuals(PRmodel1))

leveneTest(residuals(PRmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   6  1.9462 0.07911 .
##       116                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Lifestage is significant in analysis of P:R. P:R values are barely above 1 for larvae and are lowest for recruits, emphasizing the high respiratory demand not met by photosynthesis at this stage. Ratio is higher in juveniles, providing more energy for growth and calcification.

Analysis: Fused Juveniles

Respiration

Analyze respiration at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.

Build linear mixed effect model and examine for Respiration at the colony level.

Rmodel2<-lmer((R.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Rmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##              Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## Org.Number 0.037674 0.037674     1    18  6.4345 0.02067 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions are met.

qqPlot(residuals(Rmodel2))

## 38 30 
## 18 10
hist(residuals(Rmodel2))

leveneTest(residuals(Rmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  7  2.3049 0.0975 .
##       12                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Build linear mixed effect model and examine for Respiration at the polyp level.

Rmodel3<-lmer(R.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Rmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq   Mean Sq NumDF DenDF F value  Pr(>F)  
## Org.Number 0.0017085 0.0017085     1    18  5.5023 0.03065 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions met.

qqPlot(residuals(Rmodel3))

## 38 30 
## 18 10
hist(residuals(Rmodel3))

leveneTest(residuals(Rmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  1.0323 0.4579
##       12

Linear Regressions

Display linear regression of 1) colony respiration and 2) polyp respiration.

Colony respiration:

colony.R.effects <- predictorEffect(c("Org.Number"), Rmodel2) #set x axis

colony.effR<-plot(colony.R.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("purple"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-.4, 0.1))),
                   type="response", #set response scale
                   ylab=expression(bold("Colony Total Respiration")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");colony.effR

Polyp specific respiration:

polyp.R.effects <- predictorEffect(c("Org.Number"), Rmodel3) #set x axis

polyp.effR<-plot(polyp.R.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("orange"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-.4, 0.1))),
                   type="response", #set response scale
                   ylab=expression(bold("Polyp Specific Respiration")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");polyp.effR

                   #lattice=list(key.args=list(space="right",
                                              #border=FALSE, 
                                              #title=expression(bold("Temperature - Fusion")),
                                              #cex=1, 
                                              #cex.title=1)))

In fused colonies, polyp-specific respiration and total colony respiration are both significantly affected by the size of the colony (number of fused polyps). Polyp specific respiration seems to decrease in larger colonies while total colony respiration increases. This may suggest that each constituent polyp “saves” energy being part of a larger fused colony. This is very interesting to think about!

Net Photosynthesis

Analyze net photosynthesis at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.

Build linear mixed effect model and examine for Net Photosynthesis at the colony level.

Pmodel2<-lmer((P.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Pmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Org.Number 0.30307 0.30307     1    18  30.837 2.851e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions are met.

qqPlot(residuals(Pmodel2))

## 38 31 
## 18 11
hist(residuals(Pmodel2))

leveneTest(residuals(Pmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  1.1103 0.4161
##       12

Build linear mixed effect model and examine for Net Photosynthesis at the polyp level.

Pmodel3<-lmer(P.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Pmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq   Mean Sq NumDF DenDF F value  Pr(>F)  
## Org.Number 0.0054339 0.0054339     1    18  5.1573 0.03566 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions met.

qqPlot(residuals(Pmodel3))

## 27 21 
##  7  1
hist(residuals(Pmodel3))

leveneTest(residuals(Pmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  0.5616 0.7734
##       12

Linear Regressions

Display linear regression of 1) colony net photosynthesis and 2) polyp net photosynthesis.

Colony net photosynthesis:

colony.NP.effects <- predictorEffect(c("Org.Number"), Pmodel2) #set x axis

colony.effNP<-plot(colony.NP.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("purple"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-0.1, 0.8))),
                   type="response", #set response scale
                   ylab=expression(bold("Colony Total Net Photosynthesis")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");colony.effNP

Polyp net photosynthesis:

polyp.NP.effects <- predictorEffect(c("Org.Number"), Pmodel3) #set x axis

polyp.effNP<-plot(polyp.NP.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("orange"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-0.1, 0.8))),
                   type="response", #set response scale
                   ylab=expression(bold("Polyp Specific Net Photosynthesis")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");polyp.effNP

In fused colonies, net photosynthesis increases with larger fused colonies (significant), but polyp-specific net photosynthesis does not increase (not significant). This suggests that the capacity of each polyp to photosynthesize is not enhanced with fusion, but rather there is an additive effect of photosynthesis across the larger colony.

Gross Photosynthesis

Analyze gross photosynthesis at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.

Build linear mixed effect model and examine for Gross Photosynthesis at the colony level.

GPmodel2<-lmer((GP.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(GPmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Org.Number 0.55445 0.55445     1    18  28.939 4.114e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Assumptions are met.

qqPlot(residuals(GPmodel2))

## 38 36 
## 18 16
hist(residuals(GPmodel2))

leveneTest(residuals(GPmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  1.5169 0.2509
##       12

Build linear mixed effect model and examine for Gross Photosynthesis at the polyp level.

GPmodel3<-lmer(GP.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(GPmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##              Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## Org.Number 0.013236 0.013236     1    18  8.0413 0.01096 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance. Potential normality violation, revisit.

qqPlot(residuals(GPmodel3))

## 38 29 
## 18  9
hist(residuals(GPmodel3))

leveneTest(residuals(GPmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  0.9294 0.5182
##       12

Linear Regressions

Display linear regression of 1) colony respiration and 2) polyp respiration.

Colony gross photosynthesis:

colony.GP.effects <- predictorEffect(c("Org.Number"), GPmodel2) #set x axis

colony.effGP<-plot(colony.GP.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("purple"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-0.1, 0.8))),
                   type="response", #set response scale
                   ylab=expression(bold("Colony Total Gross Photosynthesis")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");colony.effGP

Polyp gross photosynthesis:

polyp.GP.effects <- predictorEffect(c("Org.Number"), GPmodel3) #set x axis

polyp.effGP<-plot(polyp.GP.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("orange"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-0.1, 0.8))),
                   type="response", #set response scale
                   ylab=expression(bold("Polyp Specific Gross Photosynthesis")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");polyp.effGP

In fused colonies, gross photosynthesis has similar results as net photosynthesis, which makes sense. Gross photosynthesis increases with larger fused colonies (significant), but polyp-specific net photosynthesis does not increase (not significant). This suggests that the capacity of each polyp to photosynthesize is not enhanced with fusion, but rather there is an additive effect of photosynthesis across the larger colony.

P:R

Build linear mixed effect model and examine for P:R ratio (gross photosynthesis : respiration). Ratios are the same when calculated for colony or polyp level, so displaying and analyzing at the polyp level here.

PRmodel2<-lmer(ratio~Org.Number + (1|Run/dpf), data=fused_oxygen)
anova(PRmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 1.0467  1.0467     1    18  0.9065 0.3537

Check assumptions of model for residual normality and variance. Assumptions met.

qqPlot(residuals(PRmodel2))

## 39 22 
## 19  2
hist(residuals(PRmodel2))

leveneTest(residuals(PRmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  7  3.9403 0.01833 *
##       12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Linear Regression

P:R ratios:

ratio.effects <- predictorEffect(c("Org.Number"), PRmodel2) #set x axis

ratio.effGP<-plot(ratio.effects,
                   lines=list(multiline=TRUE, #color lines
                              col=c("black"), 
                              lty=c(1)), 
                   confint=list(style="bands", alpha=0.3), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(-0.1, 8))),
                   type="response", #set response scale
                   ylab=expression(bold("P:R Ratio")), 
                   legend.position="top",
                   xlab=expression(bold("Number of Fused Polyps")), 
                   main="");ratio.effGP

P:R ratio does not significantly change over colony size. Juveniles at this age maintain positive P:R ratios (above 1) regardless of size. This is likely achieved by reduced polyp-specific respiration balanced by increased colony total respiration, keeping respiratory demand consistent across colony sizes.

Linear Regression Plots

Combine regression plts above into a single panel organized by metric.

fused_regressions<-plot_grid(polyp.effR, colony.effR, polyp.effNP, colony.effNP, polyp.effGP, colony.effGP, ratio.effGP, labels = c("p=0.01", "p<0.01", "p>0.05", "p=0.04", "p>0.05", "p=0.03", "p>0.05"), label_size = 16, ncol=2, nrow=4, rel_heights= c(1), rel_widths = c(1), align="h")

ggsave(filename="Mcap2020/Figures/Respiration/fused_regressions.pdf", plot=fused_regressions, dpi=300, width=10, height=16, units="in")

Conclusions

In this data, we saw a few things that are particularly interesting to think about:
(1) Respiration and photosynthesis are different between lifestages. Larval respiration and photosynthesis are low in larvae, and just maintain energetic balance (P meeting R demands).
(2) Respiration increases at recruitment (2-4 days after recruitment, 8-10 days after fertilization) and is not met by photosynthesis, highlighting energetic vulnerabilities at this stage.
(3) For older juveniles (40 days post fertilization), there is excess in photosynthesis to meet respiratory demand that is available for growth and calcification.
(4) In juveniles, those that are fused have higher respiration and photosynthesis as the size of the colony (number of polyps fused together) is larger. Photosynthesis meets respiratory demand and is in excess, providing energy for growth and calcification at this stage.
(5) Interestingly, polyp-specific respiration decreases in larger fused colonies while total respiration increases. This suggests that by joining together in a larger colony, each constituent polyp may “save” energy.
(6) Raises interesting questions for integration timing of symbiosis after settlement.